Shear strength properties of wet granular materials 
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We investigate shear strength properties of wet granular materials in the pendular state (i.e. the 
state where the liquid phase is discontinuous) as a function of water content. Sand and glass beads 
were wetted and tested in a direct shear cell and under various confining pressures. In parallel, 
we carried out three-dimensional molecular dynamics simulations by using an explicit equation 
expressing capillary force as a function of interparticle distance, water bridge volume and surface 
tension. We show that, due to the peculiar features of capillary interactions, the major influence of 
water content over the shear strength stems from the distribution of liquid bonds. This property 
results in shear strength saturation as a function of water content. We arrive at the same conclusion 
by a microscopic analysis of the shear strength. We propose a model that accounts for the capillary 
force, the granular texture and particle size polydispersity. We find fairly good agreement of the 
theoretical estimate of the shear strength with both experimental data and simulations. From 
numerical data, we analyze the connectivity and anisotropy of different classes of liquid bonds 
according to the sign and level of the normal force as well as the bond direction. We find that 
weak compressive bonds are almost isotropically distributed whereas strong compressive and tensile 
bonds have a pronounced anisotropy. The probability distribution function of normal forces is 
exponentially decreasing for strong compressive bonds, a decreasing power-law function over nearly 
one decade for weak compressive bonds and an increasing linear function in the range of tensile 
bonds. These features suggest that different bond classes do not play the same role with respect to 
the shear strength. 
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I. INTRODUCTION 

Capillary cohesion is known to influence strongly the 
strength and flow properties of granular materials. For 
example, sandcastles keep standing basically due to small 
amounts of water between sand grains [? ? ] . At low 
levels of water content, the water forms a discontinuous 
phase composed of interparticle bridges that are unevenly 
distributed in the bulk (the pendular state). Obviously, 
cohesion effects appear only at low confining pressures, 
e.g. in surface soils. It is a common observation that, 
when plowing a wet granular soil, large cohesive aggre- 
gates are formed. The largest capillary cohesion force 
for millimeter-size sand grains is about 4 x 10""' N inde- 
pendently of meniscus volume. This force is nearly four 
times the grain weight, allowing thus for the formation 
of cohesive aggregates. Transformations involving pri- 
mary particle agglomeration into coherent granules are 
of special interest in many applications in a wide range 
of industries such as pharmaceuticals, agronomic prod- 
ucts and detergents [? ? ]. 

Although capillary phenomena at the interface be- 
tween two solid bodies are well understood, it is much 
less clear how a collection of grains reacts to the pres- 
ence of a liquid. The issue is basically the same as in 
cohesionless granular media where the influence of inter- 
particle friction on the shear strength depends both on 
the features of the friction law itself and the granular 
structure. Similarly, the question here is the extent to 
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which the features of cohesion interactions are reflected 
in a global property such as shear strength. At least two 
factors seem to be important: the local force thresholds 
and the distribution of cohesive bonds in the bulk [? ? ] . 
There are reasons to think that the spatial distribution 
of water bonds should prevail over the absolute amount 
of water going to each bond. Our results, as we shall see 
in this paper, credit this point. 

There is a vast literature dealing with the experimental 
behavior of wet granular media. In soil mechanics, where 
stability considerations are of primary importance, the 
influence of water content on unsaturated earth is mostly 
studied in shear tests through the Coulomb cohesion pa- 
rameter representing the shear strength at zero confining 
stress [? ]. On the other hand, the flowability of wet pow- 
ders is expressed in terms of tensile strength measured as 
the fluidization threshold under vibration or air flow [? 
] . Direct measurements of tensile strength have been re- 
ported more recently [???]. In granular media, it is 
generally much more difficult to access local information 
such as contact forces or liquid bonds. Few investigations 
have recently been reported to visualize liquid bonds by 
means of the index matching technique [? ? ]. These 
observations underline the influence of the distribution 
of water bonds. 

A detailed description of the behavior at the parti- 
cle scale can be obtained by means of molecular dynam- 
ics (MD) simulations with an appropriate prescription of 
force laws. Recently, several simulations of wet granular 
media have been reported [???]. Mikami et al. used 
this type of simulation together with a regression expres- 
sion for the liquid bridge force as a function of liquid 
bridge volume and separation distance between particles 



[? ]. They mainly studied bubbling behavior and ag- 
glomerate formation in a fluidized bed and they found 
realistic results. Dense agglomerates were simulated by 
Groger et al. using a cohesive discrete element method 
[? ]. They found a good agreement with experimental 
data for the yield stress at all confining pressures down to 
the value of the tensile stress. Shear strength behavior of 
unsaturated granulates was also studied numerically by 
Jiang et al. as a function of suction (pressure difference 
between liquid and gas) [? ]. 

In this paper, we investigate the cohesive behavior of 
wet dense granular packings under monotonous shear- 
ing by means of experiments and three-dimensional MD 
simulations. We also analyze the shear strength from a 
microscopic expression of the stress tensor. The experi- 
ments are described in section^] We study the Coulomb 
cohesion as a function of water content for four different 
materials. The simulations are presented in section ITlTl 
We use an explicit expression of the capillary force as 
a function of interparticle distance, bridge volume and 
surface tension. The Coulomb cohesion is studied as a 
function of water content. We also investigate the bond 
connectivity and force distributions. In section IIVI we 
propose a new expression for the shear strength that ac- 
counts for particle polydispersity as well as material and 
structural parameters. We finally compare the experi- 
mental and numerical results to the theoretical predic- 
tions. 



II. EXPERIMENTS 

The experiments were designed to measure the shear 
strength at low confining pressures (< 1 kPa). The prin- 
ciples of the device are similar to those used in several 
other investigations [????]. We present here the 
setup, the materials, the wetting protocol and our main 
results. 



A. Experimental setup 

A schematic representation of the shearing setup is 
shown in Fig. ^ The wet grains are poured in a plexiglas 
cylindrical cell and confined by means of a circular lid of 
area S placed on top of the material. The lid is equipped 
with a reservoir allowing to impose an overload by adding 
desired amount of sand. The total vertical force A'' acting 
on the sample is the sum of the weights of lid and sand 
(shown by A on the figure) . The cell is composed of two 
disjoint parts kept together during sample preparation. 
The upper part can move horizontally with respect to the 
lower part by pulling on a rope attached to it and which 
supports a cupel through a pulley (shown by B on the 
figure). The pulling force T can be increased by adding 
sand into the cupel. The friction force between the two 
parts of the cell is reduced by water lubricating the rims 
and we checked that it remains negligibly small during 
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FIG. 1: Testing cell and shearing setup. 



shear. In order to reduce the friction force exerted by 
the material along the walls (Janssen effect), the thick- 
ness h of the upper part of the sample was taken to be 
below the diameter of the cell (46 mm). The heights of 
the upper and lower parts are about 10 mm and 15 mm 
respectively. The sample is sheared along the common 
section of the two parts of the cell. This shear plane is 
subjected to a tangential stress t — T / S and a normal 
stress (T = N / S + pgh, where p is the bulk density and g 
is the gravity. 

In the experiments, we gradually increase the shear 
stress T for a fixed value of a. Unstable shearing occurs 
when T reaches the shear strength r,„ , resulting in a sud- 
den slide of the upper part of the sample. The upper part 
is stopped by collision with two bars located 5 mm away 
from the cell. We did not measure the displacements. We 
recorded t„j for different values of a in the range varying 
from 200 Pa to 800 Pa, and for different values of water 
content. 



B. Materials and wetting protocol 

We used four types of materials: (1) a sand composed 
of angular grains with diameter in the range from 0.1 mm 
to 0.4 mm, (2) "tightly-graded" polydisperse glass beads 
with diameters from 0.4 mm to 0.5 mm, (3) "well-graded" 
polydisperse glass beads with diameters from 0.4 mm to 
0.8 mm, and (4) monodisperse glass beads of diameter 
1 mm. 

In order to wet the grains, we add distilled water to 
dry material placed in a vessel which is then closed and 
energetically shaken for about one minute. The vessel 
used for mixing is transparent so that during shaking we 
can check visually whether the water is homogeneously 



mixed with the grains. In particular, we continue shaking 
until all visible water clusters disappear. Increasing the 
duration of shaking beyond 1 minute did not change the 
measured values of the shear strength. After mixing, the 
wetted material is poured into the testing cell. The water 
content is evaluated by comparing the masses of a sample 
of the material before and after testing by means of a 
heat chamber used for drying the sample at 105°C. The 
water content is given by w = mwjms^ where rriw and 
mg are the masses of water and grains, respectively. The 
wet materials were tested for water contents below 0.05 
corresponding to the pendular state for our materials. 
The experiments were performed at ambient conditions. 
Each experiment lasted a few minutes. The loss of liquid 
was always below two percent. This loss is not only due to 
evaporation but also due to partial wetting of the internal 
walls of the cell. It is small enough to assume a constant 
liquid volume (as in simulations, see below). 



C. Results 

Figure I^Ja) shows the yield loci t-ct for the sand at 
several levels of water content w. Within experimental 
precision, the data are well fitted by a straight line, in 
agreement with the Mohr- Coulomb model 



T = /i(T 



(1) 



where /i = tan ip is the internal coefficient of friction and 
c is the Coulomb cohesion. It is remarkable that ip ~ 
33° is almost independant of w. On the other hand, 
the Coulomb cohesion increases nonlinearly with w and 
saturates to c,„ ~ 600 Pa at w™ ~ 0.03, as shown in 
Fig.Hb). 

For tightly-graded polydisperse beads we get a simi- 
lar behavior with <p ~ 30°; Figs. ISJ^a) and (b). How- 
ever, saturation occurs at a lower level of water content 
(wm — 0.025) and cohesion (c„i ~ 350 Pa). Figures 
2[a) and (b) show the results for well-graded polydis- 
perse beads. Saturation occurs at about Wm — 0.01 
and Cm — 300 Pa. Enhance fluctuations observed in 
the data from glass beads, compared to the sand, may 
be attributed to a lower level of cohesion and a tighter 
size distribution of glass beads. In fact, a lower level 
of cohesion leads to larger mobility of the particles and 
the porosity increases as the particle sizes are less widely 
distributed. In the case of monodisperse beads, the co- 
hesion jumps from zero for the dry material to a nonzero 
value {cm — 150 Pa) independently from water content; 
Figs.inia) and (b). Since we have few data points between 
w ^ and w ~ 0.01, the saturation level Wm should be 
below 0.01. 

Notice that the shear tests provide quite reproducible 
results. We see that, in Figs. |2Ia)-[3Ia), the Mohr- 
Coulomb line passes through most of data points. The 
very weak dispersion of the data points about this line 
shows the high reproducibility of the testing procedure. 
This means that the experimental measurement of the 
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FIG. 2: (a) Yield loci r-a of sand for increasing level of water 
content; (b) the Coulomb cohesion as a function of water 
content. The dashed line is drawn as a guide to the eyes. 



Coulomb cohesion is reliable. The different values of c^ 
for different materials will be discussed in section llVl fsee 
also Table Pi. 

The value of Wm is less clearly defined and it is likely 
to depend on two factors: (1) the surface state of the 
particles and (2) possible clustering of the liquid phase 
[???]. The sand grains have a rough surface requiring 
more water to form a meniscus than glass beads which 
are much more smooth. On the other hand, partial clus- 
tering of water may occur and this might require a larger 
amount of water for the formation of liquid bridges al- 
though we observed no clustering at the visible parts of 
the packing through the transparent walls of the testing 
cell. It is worth noting that it is not straightforward to 
evaluate the variability of w since each value results from 
several experiments. If the evaluation is based only on 
the measurement of the water content at the beginning 
and at the end of each test, the error would be as small 
as 2% and this can not be shown in the figures (it would 
be smaller than the size of the data point symbols). 
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FIG. 3: (a) Yield loci r-a of tightly-graded glass beads for 
increasing level of water content; (b) the Coulomb cohesion 
as a function of water content. 



FIG. 4: (a) Yield loci r-a of well-graded glass beads for in- 
creasing level of water content; (b) the Coulomb cohesion as 
a function of water content. 



III. SIMULATIONS 

For the simulations, we employed the framework of 
the molecular dynamics method [? ? ]. This method 
is referred to as Distinct Element Method (DEM) in the 
geotechnical context. The heart of our simulations is, 
however, the model of capillary cohesion and its imple- 
mentation, as well as the way liquid bridges are numer- 
ically distributed in the packing. These aspects are de- 
tailed below. 



act also through a Coulomb friction law with a viscous 
regularization at low sliding velocities. The material is 
sheared quasistatically in a direct shear cell in the pres- 
ence of gravity and confining stresses. The damping pa- 
rameter and the normal stiffness (force per unit overlap) 
were adjusted in order to get largest time step (10~^ s) 
and small overlaps within numerical stability. The nor- 
mal stiffness and the interparticle coefficient of friction 
were 10"^ N/m and 0.4, respectively. 



A. Molecular Dynamics 

We implemented the basic molecular dynamics method 
for spherical particles. The equations of motion are inte- 
grated according to the velocity Verlet scheme. The nor- 
mal force between two particles is the sum of a repulsive 
force as a linear function of the overlap and an attraction 
force due to the presence of a liquid bridge at contact or 
for a gap up to a rupture distance (see below). As usual 
in molecular dynamics simulations, normal dissipation is 
accounted for by viscous damping. The particles inter- 



B. Capillary cohesion law 

The capillary attraction force between two particles 
is a consequence of the liquid surface tension and the 
pressure difference between liquid and gas phases [? ]. 
For efficient numerical calculation, we need an explicit 
expression of the capillary force /^ as a function of the 
interparticle gap J„. By extending the work of Mikami 
et al. [? ], it was recently shown by SouHe et al. [? ] 
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FIG. 5: (a) Yield loci r-cr of monodisperse glass beads for 
increasing level of water content; (b) the Coulomb cohesion 
as a function of water content. 



that the capillary force can be cast in the following form: 
-7^7,^/i^^R;{exp(A5;^ + B) + C} for (5* > 

for St < 



/, 



-Tr-f,^Wlh{exp{B) + C} 



(2) 

where i?i and i?2 are the sphere radii (i?i < R2), 7s is 
the liquid surface tension, J* = Sn/R2 (see Fig. El. The 
parameters A, B and C are functions of the liquid volume 
Vb of the bond and the contact angle as follows 

A = -l.l(Vb*)-"-53, 



B = (-0.1481nT4*-0.' 



where V^^ 



C -- 
Vb/Rl 



-0.0082 lny* + 0.48, 



0.00181nK* + 0.078, 



(3) 



It can be shown that a liquid bond is stable as long as 
the gap is below a de-bonding distance 5™°^ given by [? 



(1 -f 0.561) Vb^/^ 




FIG. 6: Typical behavior of the capillary force /^ as a function 
of the gap Sn (solid line). The elastic repulsive normal force 
as a result of overlapping {&„ < 0) is shown (short-dashed 
line). The resultant normal force in this range is shown as 
well (long-dashed line), z'^'^ , z'^~ and z~ are the partial wet 
coordination numbers in this range. Inset: Geometry of a 
liquid bridge between two particles. 



Figure El shows a schematic representation of the capil- 
lary force as a function of the gap. In Fig.[7|the capillary 
force /^ is displayed as a function of the gap (5„ accord- 
ing to Eq. 121 for different values of the liquid bridge vol- 
ume Vfe. These plots are in perfect agreement with those 
obtained by other authors by integration of the Laplace- 
Young equation and verified experimentally [? ]. The 
largest absolute value /o of the capillary force occurs at 
(5„ = 0. It is remarkable that /o is directly proportional 
to y/RiRa and only very weakly dependent on the liq- 
uid volume. This property, which might seem counter- 
intuitive, is important for the model that will be intro- 
duced in section Hvl Hence, with a good approximation, 
we may write 



/o = Ky/RlR2, 



(5) 



(4) 



where k is a function only of the surface tension and the 
contact angle. In our case, with glass beads and water 
bridges, we have k = 0.4 N/m. Let us underline here the 
fact that the capillary bond in the range Sn > is unsta- 
ble with respect to the forces acting on two particles. In 
other words, when pulling two particles apart from one 
another, the liquid bond fails at zero gap for /„ = — /o- 
In our simulations, we find that the fraction of liquid 
bonds in the range Sn > is always below 15% (see Fig. 
Illfl . This shows that the capillary failure threshold /o 
is far more important for the failure of a wet material 
than the de-bonding distance (5™°^. This point will be 
discussed in more detail in section Hvl 
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FIG. 7: Capillary force as a function of the gap according to 
Eq. 121 for increasing liquid volume Vj,. 



C. Sample preparation 

The numerical samples are com.posed of spherical par- 
ticles of three different diameters (2 mm, 1.5 mm and 1 
mm) placed randomly in a cylindrical cell in appropri- 
ate proportions (50%, 30% and 20%) to represent one of 
our experimental samples composed of glass beads. The 
initial configuration is prepared under gravity without 
introducing capillary bonds. Then, we attribute a capil- 
lary bond to each pair of particles within the de-bonding 
distance. Finally, the sample is consolidated under the 
action of a vertical confining pressure with a zero coef- 
ficient of friction. The consolidation is stopped and the 
coefficient of friction set to 0.4 as soon as the solid frac- 
tion (\) = 0.6 is reached. The subsequent compaction is 
negligibly small. 

The volume Vf, attributed to a capillary bond between 
two particles is taken to be proportional to the particle 
diameters and the intercenter distance, and such that the 
total volume of all liquid bonds in the sample is equal to 
that of the added water. Since all particle pairs within 
the de-bonding distance are considered, the liquid coordi- 
nation number z (i.e. the average number of liquid bonds 
per particle) obtained by this procedure has the highest 
possible value. For our sample we get z — % for w — 0.01 
(see Fig. ll4|l . In our simulations, the liquid bond volumes 
vary by a factor 8 from the contact between the small- 
est particles to that between the largest particles. On the 
other hand, the contact angle with a good approximation 
was set to zero. Moreover, since the de-bonding length 

1/3 

varies as Vj, , there is a factor 2 between the shortest 
and longest de-bonding distances. 

During shearing, the number of liquid bonds evolves 
and the available liquid must be redistributed in the sys- 
tem. We used two different methods for redistribution: 
(1) we simply apply the above procedure every time the 
contact list is updated; (2) the volume of a broken liquid 
bond is split between the corresponding particles (pro- 



FIG. 8: The shear stress r as a function of shearing distance 51 
normalized by the average particle diameter (d) for a dry and 
a wet sample simulated by the molecular dynamics method. 
The confining stress is a = 300 Pa. 



portionally to their diameters) and conserved for the for- 
mation of new liquid bonds when a contact occurs with 
the same particles. In this method, the volume of free 
liquid left after de-bonding is kept with the two particles 
(and not distributed to the other bonds of the same par- 
ticles) and used only if a new contact is formed. This 
implies that, if the initial liquid distribution is homoge- 
neous, then it will remain so during deformation as in 
the first method. In other words, the liquid will not mi- 
grate considerably and one should expect quite similar 
results from both methods. Indeed, in different tests, we 
found that both methods lead to nearly identical results. 
Unless mentioned explicitly, all results presented in this 
section were obtained by the first method. 



D. Boundary conditions and driving 

As in experiments, the cylindrical cell is composed of 
two disjoint parts. The lower part is fixed whereas the 
upper part moves horizontally, giving rise to a shear plane 
along the common section of the two parts. We apply a 
constant vertical load a, the same as in experiments, on 
top of the sample. However, in contrast to experiments, 
shearing is controlled by imposing a constant horizontal 
velocity on the upper part. The numerical sample has 
exactly the same dimensions as in experiments. 



E. Results 

Figure |21 shows the stress-strain plot for a dry and a 
wet sample with w = 0.01. The initial configuration is 
the same in both simulations. The initial elastic increase 
of r (up to ^ 75 Pa) as a function of 5(, is common be- 
tween the two samples. We observe no stress peak in the 
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FIG. 9: (a) Yield loci r-cr from 15 simulations for increasing 
level of water content; (b) the corresponding values of the 
Coulomb cohesion c as a function of water content. 



FIG. 10: (Color online) Compressive force network (a) and 
tensile force network (b) in a thin vertical layer. Weak forces 
are not shown, xz is the shear plane and j/-axis points upward. 



wet case. The steady state ("critical state" in soils me- 
chanics) is reached at Si ~ (d) for all water contents. The 
steady state deformation involves numerous instabilities 
that occur throughout the system and appear in the form 
of rapid stress drops on the stress-strain plots. We see 
that in transition from dry to wet materials the frequency 
of such instabilities declines while their amplitudes grow. 

We now turn to the evolution of the Coulomb cohesion 
as a function of w. Figure El^a) shows fitted yield loci 
from 15 simulations involving three different values of 
the confining pressure a and five different values of the 
water content w. The Coulomb cohesion c is drawn as 
a function of w in Fig. El^b). The latter is very similar 
to the corresponding experimental plot (Fig. EJb)) for 
monodisperse glass beads. We observe a saturation of 
c at still lower levels of water content {wm — 0.001). 
Note also that, while the average grain size is nearly the 
same in these simulations and in the case of monodisperse 
experimental glass beads, the maximum cohesion c^ = 
120 Pa in the simulations is very close to that (150 Pa) 
for 1 mm glass beads. 

In contrast to the experiments, where the stresses are 
measured at the walls, it is also possible to compute the 
stress tensor for grain-to-grain forces in the simulations 
(see below). However, we found that in our simulations, 
the results from these two methods do not coincide. This 
is because in direct shearing, wall effects give rise to large 



stress gradients in the bulk. This point has been ana- 
lyzed in detail by Thornton and Zhang [? ]. We used 
here the wall stresses r and a for comparison with the 
corresponding experimental values. 

Figure ^| displays a typical example of the force net- 
work in a thin vertical layer (parallel to xy plane). We 
observe a strongly inhomogeneous transmission of both 
compressive and tensile forces. The liquid bonds belong 
to three different classes: (1) contacts carrying a com- 
pressive (positive) force; (2) contacts carrying a tensile 
(negative) force; (3) liquid bonds with no contact and 
thus carrying a tensile force. We denote the correspond- 
ing partial coordination numbers by z*^"*", z'^~ and z~, 
respectively (see Fig. ^ . Fig. ^2 shows the evolution of 
these partial coordination numbers as a function of shear 
displacement. Interestingly, although the bonds are op- 
timally distributed in the sample (according to the first 
redistribution method), z~ is only about one bond per 
particle. This means that the overall contribution of this 
class, involving a low number of bonds and rather weak 
forces, to stress transmission is marginal. The evolution 
of particle connectivity is mainly refiected in the regular 
fall-off of z^~ during shear. This class, with numerous 
contacts and a high force level, provides the largest con- 
tribution to the transmission of tensile stresses in the 
packing. The rather large and constant value of z'^'^ is 
a reflection of the fact that the packing is globally sub- 
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FIG. 11: Partial coordination numbers as a function of shear- 
ing distance. 



jected to boundary compressive stresses. 

The probability distribution function P of normal 
bond forces is shown in Fig. E| The largest forces be- 
long to the compressive network (involving no limiting 
threshold) whereas the tensile forces extend down to the 
capillary force threshold — /q. The distribution reveals 
three different force intervals with different statistics in 
each interval: 



P{fn) « 



fn±f0 I 



for /„ > (/«) 
for < /„ < (/„) 
^0 for - /o < /„ < 



(6) 



where (/„) is the average normal force, q;~0.4,/3~0.16, 
7 ~ 0.15 and Pq = P{—fo). The distribution for com- 
pressive forces is reminiscent of that observed in dry gran- 
ular media [???]. The exponent P has, however, a 
smaller value in the cohesive case. This means that ten- 
sile forces allow the packing to sustain stronger compres- 
sive force chains than in a cohesionless material. The 
power law distribution of weak compressive forces (the 
range < /„ < (/«)) over nearly one decade suggests 
that, in the presence of cohesive bonds, the contact net- 
work in this range is self-similar. These contacts are re- 
ferred to as "weak contacts" and they play an important 
role in propping strong force chains in the complemen- 
tary "strong network" (the range /„ > (/«)) [? ]. The 
probability distribution P is robust and the exponents 
do not evolve during shear. 

Figure [T^ a^ shows a polar diagram of the probability 
distribution function P{n) of bond directions n at the 
end of shearing for w — 0.02 (the distribution being sim- 
ilar in other cases). The distribution is nearly isotropic 
along the shear plane (xz plane in the figure) but it shows 
a pronounced anisotropy in the shear direction along the 
vertical plane (xy plane in the figure). Fig. ll^f bl displays 
three separate polar diagrams of the strong, weak and 
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FIG. 12: Probalility distribution function P of normal forces 
/„ normalized by the mean (/„) represented on linear scale 
(a), log-linear scale (b) and log- log scale (c). 



tensile bonds along the vertical plane. The strong com- 
pressive bonds form an anisotropic distribution with its 
longest axis oriented at 45° to the horizontal, as expected. 
The weak compressive bonds have a nearly isotropic dis- 
tribution. Finally, the tensile bonds have an anisotropic 
distribution whose longest axis is perpendicular to that 
of strong compressive bonds. It has also been argued 
that the tensile forces play the same role in sustaining 
the strong force chains as the weak network [? ]. 




FIG. 13: (Color online) (a) Polar diagram of the probability 
distribution function of bond direction; (b) polar diagrams of 
strong (solid line), weak (dashed line) and tensile (dotted line) 
bonds along the vertical plane. The direction of extension and 
compression are represented (arrows). 



IV. A MICROSCOPIC ANALYSIS 

In this section, we would like to propose an expression 
for the Coulomb cohesion c as a function of parameters 
pertaining to the granular microstructure in the presence 
of Hquid bonds. When the Mohr-Coulomb model (Eq.P) 
is valid also in the range of negative stresses down to the 
tensile strength — cr*, the Coulomb cohesion is related 
to the tensile strength by c = /icr*. Keeping with this 
assumption, we therefore consider the tensile strength. 
Theoretical evaluation of the tensile strength and its com- 
parison with experiments or simulations can be found in 
recent literature concerning wet granular materials [? ? 
] . The first theoretical expression of tensile strength was 
proposed by Rumpf [? ] . The key point is how capillary 
forces are mobilized in a wet material and what are the 
relevant structural parameters. Several parameters, such 



as the particle size and its distribution, the solid fraction 
(/) and the internal coefficient of friction fj,, influence the 
macroscopic cohesion. The effect of the water content w 
is somehow counter-intuitive. Indeed, the Coulomb co- 
hesion saturates as the water content is increased while 
we expect that both the number of liquid bonds and the 
strength of each bond should increase with the water con- 
tent and lead to higher cohesion. 

In order to estimate the tensile strength from contact 
forces, we consider the general expression of the stress 
tensor cr in a granular material. This is an average quan- 
tity with a well-established expression involving contact 
forces /'"' and inter-center distances {.''[? 7 ]: 






(7) 



kev 



where i and j refer to components and V is the control 
volume. The derivation of this expression is independent 
of the nature of interactions so that Eq.[7|holds also in the 
presence of capillary forces. In this case, the set of con- 
tact points is simply extended to cover capillary bridges. 
It can be shown that cr is symmetric if no torques are 
transmitted at the interaction points. From Eq. \7\ the 
stress CTii in the direction of extension is given by: 
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(8) 



where Uw is the number of bonds per unit volume, /i and 
£i are the components along the direction of extension. 
The symbol (• • •)j designs averaging over all bonds in 
the control volume V. Let n and t be the normal unit 
vector and a tangential unit vector, respectively, at a 
given bond. Then, /i = /„ni + ftti and ii = £ni where 
/„ and ft are the normal and tangential components of 
the force, and £ is the length of the branch vector £. Let 
us set rii = cos 9 and ti — sin 9, where 9 is the angle 
between n and the direction of extension. Substituting 
in Eq. |S1 and assuming for simplicity that /„ and ft are 
independent of 9, we arrive at the following expression: 



CTll = - n^ifn i)b- 



(9) 



In solid state physics, a theoretical tensile strength cr"' 
is introduced from inter-atomic forces by assuming that 
the same failure threshold is reached simultaneously for 
all pairs of atoms in the direction of traction [? ]. In a 
similar approach, we may introduce a theoretical tensile 
strength for a wet particle assembly by replacing /„ in 
Eq. Oby the capillary force threshold /q: 
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c^"" = 2 '^'"^■^° ^^'>- 



(10) 



The bond density n.w is simply half the average number 
of bonds per particle divided by the free volume, i.e. the 
mean volume 1^ of a Voronoi cell surrounding the par- 
ticle. The latter is simply the average particle volume 
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(l/6)7r((i^) divided by the solid fraction (j). Introducing 
these expressions in Eq. IIUI and using Eq. |S1 we get 
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where z is the average number of bonds per particle and 
we have 



(dl/2)(d)(d3/2 
W> 



(12) 



In derivation of the expression of s, it was assumed that 
the particle radii i?i and i?2 are not correlated. This 
means that various granulometric classes are homoge- 
neously distributed in the bulk. It is easy to see that 
for a uniform size distribution, s varies from 8/15 to 1 
as the smallest particle size increases from to the mean 
particle size. For a monodisperse assembly, we have s = 1 
(see TablePl. Eq.^Jis similar to the expression proposed 
first by Rumpf [? ] for monodisperse materials (so, with- 
out the s prefactor) , and recently derived from the stress 
tensor by Groger et al. [? ]. Our equation 1 1 1 1 accounts in 
a simple way for polydispersity and the correlation be- 
tween the capillary force threshold /o and the particle 
size d. For real polydisperse materials, the factor s is 
crucial for comparing the model with experiments. 

By definition, the Coulomb cohesion is the yield shear 
stress at zero confining pressure in which case the cap- 
illary forces are the only forces acting in the material. 
In this limit, the normal stress a on the shear plane is 
simply equal to the average capillary force divided by 
the sample section S. We have seen that "gap" liquid 
bonds (without contact) contribute only marginally to 
force transmission (z~ being below one bond per parti- 
cle); see Fig. ^2 It is thus reasonable to assume that 
the capillary force at each bond is /q. This means that, 
in the absence of confining stresses, we may set a = a*^. 
Then, the shear stress at yield is the theoretical cohesion 
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It is important to note that the water content does not 
enter the above expression of c"*. The only parameter 
related to water is k. This suggests that the water con- 
tent manifests itself mainly through the wet coordination 
number z = z^ + z'^~ + z'^'^ . In particular, the cohesion 
Cm at saturation corresponds to the saturation of z as 
the water content w is increased. In fact, when a cer- 
tain amount of water is homogeneously distributed in the 
whole sample within the de-bonding distance, one finds 
that z increases with w and saturates beyond w = Wm- 
This is shown in Fig. E| for the initial configuration of 
our numerical samples. This is a purely geometrical ef- 
fect related to steric exclusions among particles and it 
explains therefore the saturation of cohesion with water 
content according to Eq. E| 
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FIG. 14: The initial wet coordination number in simulations 
as a function of water content. 



Although this geometrical saturation should dominate 
in the pendular state, it does not elude that other mecha- 
nisms might play a role in the experiments. In particular, 
as the liquid content is increased, the liquid bonds may 
coalesce at least locally, leading to bond saturation. We 
did not observe liquid bond clustering at the visible parts 
of our samples. 

The initial value of z depends also on the prepara- 
tion protocol. Since /o is independent of w, the same 
amount of water can be distributed in such a way as to 
produce a lower number of liquid bonds and thus a lower 
macroscopic cohesion. This effect is illustrated in Fig. [Tsl 
where the stress-strain plots are shown for two initially 
identical configurations differing only in the number of 
liquid bonds for the same water content. The simulation 
was carried out by using the second method of liquid re- 
distribution (section lllll) . In the sample where half of 
the water bonds has been removed, the wet coordination 
number increases with deformation. But in the initial 
stages of deformation, the cohesion is close to half that 
of the sample involving a double number of water bonds, 
and it increases as the wet coordination number grows. 

We summarize in Tabled the theoretical estimates c*'' 
and the measured values c,„ of the saturated Coulomb 
cohesion for all our experimental and numerical samples 
together with the values of the parameters involved in 
Eq. El The value of the polydispersity factor s was cal- 
culated from the knowledge of the particle size distribu- 
tions. Note that c*'* is in excellent agreement with Cm 
both for experiments and simulations. It is noteworthy 
that if the prefactor s were not incorporated in Eq. [T51 
(i.e. if Rumpf 's equation had been used), the measured 
value of Cm would be below c*'' by a factor s. 

This agreement between the theoretical estimate and 
experiments with sand and glass beads was obtained with 
z — 6 which is a reasonable value of the bond coordina- 
tion number in the case of a homogeneous distribution of 
water in the bulk, and it is suggested by recent experi- 
mental observations [? ? ]. Note, however, that in the 
case of fine- grain samples (sand and GBl in Table QJ a 
closer agreement can be obtained with a lower value of z. 
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FIG. 15: The shear stress r as a function of shearing distance 
(5£ normalized by the average particle diameter (d) for a dry 
and two wet samples with different numbers of liquid bonds; 
see text. The confining stress is ct = 300 Pa. 



TABLE L Measured and theoretical parameters for all our 
experimental and numerical samples. The approximate value 
of the bond coordination number z for the experiments (indi- 
cated by a question mark) was suggested by the literature [? 





Sand 


GBr 


032' 


GB3'= 


Simulations 


(d> (mm) 


0.16 


0.45 


0.60 


1.00 


1.65 


s 


0.50 


0.99 


0.91 


1.00 


0.79 


z 


6(?) 


6(?) 


6(?) 


6(?) 


9 





0.6 


0.6 


0.6 


0.6 


0.6 


M 


0.66 


0.58 


0.58 


0.46 


0.48 


Om (Pa) 


600 


350 


300 


150 


120 


c"' (Pa) 


709 


438 


302 


158 


118 



"tigthly-graded polydisperse glass beads 
''well-graded polydisperse glass beads 
'^monodisperse glass beads 



This is suggestive in the sense that liquid bond clustering 
might indeed occur more frequently for fine grains and 
reduce thus the effective bond coordination number. 

By construction, the theoretical estimate is an upper 
bound for cohesion. For brittle materials, failure is ini- 
tiated by the breakdown of a few bonds and propagates 
subsequently into the material. When this mechanism 
works, the tensile strength is not controlled by the aver- 
age stress (as assumed in the derivation of ct*'*) but by the 
largest local stresses, and hence the effective strength is 
far below the theoretical one (in proportion to the stress 
concentration factor) [? ]. Hence, the nice agreement 
of the theoretical estimate both with simulations and ex- 
periments suggests that the failure of our wetted granular 
materials is ductile and the shear strength is controlled 
by the mean tensile force. 



V. CONCLUSION 

We performed experiments and discrete element sim- 
ulations to analyze the Coulomb cohesion of wet granu- 
lar media in the pendular state. It was shown that the 
Coulomb cohesion increases with water content and sat- 
urates to a maximum value that depends only on the 
nature of the material. An interesting aspect that was 
partly investigated by simulations is that the cohesion 
is basically controlled by the number of liquid bonds. 
This suggests that the saturation of Coulomb cohesion 
occurs since new bonds are hardly formed beyond a cer- 
tain amount of the water content. On the other hand, 
the capillary failure threshold is nearly independent of 
the local liquid volume. 

Starting with the expression of the stress tensor, we 
also introduced a novel expression for the Coulomb cohe- 
sion as a function of material and structural parameters. 
This expression extends the classical model of Rumpf to 
polydisperse materials. We found that our model is in ex- 
cellent agreement with experimental and numerical data. 

From numerical data, we analyzed the connectivity and 
anisotropy of different classes of liquid bonds according 
to the sign and level of the normal force as well as the 
bond direction. We found that weak compressive bonds 
are almost isotropically distributed whereas strong com- 
pressive and tensile bonds have a pronounced anisotropy. 
It was shown that the probability distribution function 
of normal forces is exponentially decreasing for strong 
compressive bonds, a decreasing power-law function over 
nearly one decade for weak compressive bonds and an 
increasing linear function in the range of tensile bonds. 

In the extension of this work, it is essential to evaluate 
the limits of the model by considering other materials 
and non monotonous loading paths. In particular, we 
would like to study the shear strength of granular me- 
dia with a larger polydispersity than materials that were 
used in the present investigation. Since the distribution 
of liquid bonds seems to be a major parameter for the 
cohesion, it also merits to be investigated in more detail 
experimentally. Finally, an interesting application of the 
ideas put forward in this paper would be to examine by 
which mechanisms the cohesion of a sample of wet sand 
increases as a result of compactification. 
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